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1.  Introduction 

Tl'.e  advent  of  large  scale  digital  computers,  with  their 
phenomenal  speed  and  accuracy  of  computation,  has  created  a 
challenge  to  mathematicians  to  use  the  full  capalllltles  of 
the  machinery  to  advance  knowledge  wltnln  their  field.  It  Is 
In  this  spirit  that  this  paper  has  teen  written. 

Nonlinear  dynamics  proMems  prove  difficult  to  Investigate 
analytically  or  numerically  If  probabilistic  terms  are  Involved 
In  their  equations.  This  paper  develops  a  technique  for  hand¬ 
ling  problems  of  this  type  that  Is  practical  only  If  a  high 
speed  computer  Is  avallatle  to  perform  the  attendant  compu¬ 
tations  . 

The  method  Is  general,  lut  will  be  applied  to  the  Van  der 
Pol  equation: 

X  -♦>  X(x^  -  1  )x  +  X  -  r(  t ) , 

(1) 

x-u,  x-v,  at  t-t^. 


*The  RAND  Corporation. 
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(A  dotting  indicate!  differentiation  with  respect  to  t.) 

The  classical  problem  considers  r(t)  ■  A  cos  wt.  We 
consider  r(t)  to  be  a  random  function  with  a  known  probaba- 
billty  distribution. 

2.  Theory 

We  first  use  the  phase  plane  form  of  the  equation: 


X  -  y. 

(2) 

y  ■  —  X(x^  -  l)y  -  X  +  r(t). 

Let  x(t^)  -  x^,  y(t^)  -  y^,  r(t^)  -  r^,  ^ 

where  x(tQ)  ■  u,  y(tQ)  ■  v.  The  Euler  solution  procedure 
for  the  foregoing  equation  Is: 


(3) 


(4) 


I 


*n+l  ■  ='n  K 


^n+l  ■  ’‘(*0  -  "  ’'n 


Let  us  now  define  a  probability  function 


P,  -  P(C(xN,  y^);u,v). 


the  probability  that  at  time  t^  the  Integral  curve  with 
initial  conditions  Xq  -  u,  yQ  “  will  satisfy  the 
condition  0(x^»  ^n^’ 

It  is  assumed  that  Pq  is  known  at  all  points  of  the  x,y 
phase  plane.  Then  if  r  Is  a  fixed  variable, 
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’’n+l  '  Wl): 


■  K-  ^(’'n  "  "  ’'n  ^ 


P(C(Xn,  y^(r));u,v). 


Since  r(t)  has  a  given  prol  ability  distribution,  0(r), 
we  must  average  the  right-hand  side  over  r.  This  may  be 
expressed  as 


(5) 


'’n+l^Vr  yn+l  =  “''')  ■c/’f<’'n-  Yp  >-)  • 


Let  this  fonrallsm  now  be  applied  to  a  question  of 
stability.  This  requires  an  Investigation  of  for  large  n, 

where 


Mx^I  <  a, 

j 

[lypl  <  b- 


If  a  and  b  are  chosen  sufficiently  large,  P^  will 
converge  to  1  uniformly  throughout  the  plane  provided  that  the 
solution  Is  stable. 

If  a  and  b  are  varied  for  different  runs  to  determine 

those  values  of  a,  I  at  whlcbi  P  first  becomes  less  tb.an  1 

n 

for  large  n,  a  region  may  be  set  up  In  whlcli  tbie  Integral 
curve  must  lie.  This  may  te  refined  to  an  outer-tangentlal 
rectangle  by  considering  C(x^,y^)  to  be 


r 


P-.122b 


It  may  be  suspected,  or  desired,  that  the  Integral  curves 
lie  within  a  given  region.  To  test  this,  one  assumes  the 
condition  C(x,y)  t<^  be  f(x,y)  <  0,  e.g., 

^(x,y)  -  X  y  -  R  ,  Circular  region, 

2  2?  2 

f(x,y)  -  (x  -  A  )(y  —  B  ),  Rectangular  region, 

f(3t,y)  -  (x^  +  y^  -  R^)(x^  y^  -  Annulus. 

Again,  one  tests  for  converging  to  1  for  large  n, 

Xinlformly  over  the  plane. 

In  practice,  the  Infinite  (x,y)  plane  Is  replaced  by  a 
finite  cross  grid.  Suitable  adjustments  must  be  made  to  define 
the  value  of  the  Integral  In  the  neighborhood  of  the  boundary. 
'Rie  computational  and  theoretical  errors  associated  with  this 
procedure  are : 

1.  Truncation  error  of  the  Euler  process:  This  can  be¬ 
come  very  bad  In  the  relaxation  region  and  transition  regions, 
for  large  A. 

2.  Computational  Roundoff:  Hils  becomes  larger  for 
smaller  values  of  A,  tut  Is  not  a  function  of  grid  size. 

5.  Intragrid  Interpolation  effects  and  grid— edge  errors. 

4.  Truncation  error  associated  with  the  numerical  Inte¬ 
gration  . 

Neither  the  effects  of  these  errors,  nor  their  hounds, 
have  been  Investigated  In  this  paper. 
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IlluatratlonB 

To  Illustrate  the  method  described  In  this  paper,  equation 
(l)  was  considered,  with  r(t)  assuming  the  values  ^  k  and 
—  k  with  equal  probability.  chosen  as 

Xn  <  a,  <  b. 

Two  values  of  X  were  considered,  X^-fl,  x*  —  !. 

It  was  expected  that  the  nature  of  the  results  for  small  k 

would  not  vary  too  much  from  k  -  0. 

Por  k  -  0,  one  has  the  standard  Van  der  Pol  equation.  For 

X  ■  1,  there  exists  one  periodic  solution  towards  which  all 

solutions  converge  rapldly^^^.  The  graph  of  this  periodic 

(2) 

solution  Is  shown  In  Figure  1'  .  For  a  •  b  •  4,  P^  should 

n 

converge  to  1  for  all  (x,y). 

The  value  x  -  -  1  was  selected  to  exemplify  a  more  varied 
probability  distribution.  If  In  equation  (l),  the  substitution 
■  -  t  Is  made,  equation  (l)  as  a  function  of  T  would  become 

X  -  X(x^  -  l)x  X  -  r(f). 

Por  negative  X,  this  equation  Is  Identical  to  the  previous 
case.  Hence,  the  same  unique  periodic  solution  exists,  but  with 
the  y  orientation  reversed.  However,  as  a  periodic  solution  to 
the  t  equation.  It  Is  unstable.  The  singular  point  at  the 
origin  becomes  stable. 


(1) .  Brock,  P.,  "Methods  In  Nonlinear  Vibrations,"  M.  S. 

Thesis,  N.Y.U.,  1947. 

(2) .  Andronow,  A.  A.,  Chaikin,  C.  B. ,  "Theory  of  Oscillations," 
Princeton,  1949,  p.  251. 
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Figure  2  l8  a  plot 

(x®  -  l)y  -  X  -  0. 

These  are  Isoclines  of  horizontal  slope  for  the  fcunlly  of 
Integral  curves.  The  direction  of  slopes  In  regions  A  —  P 
are  shown  in  the  figure.  For  (a,b)  chosen  as  (l,l)  and 
large  n.  It  is  clear  that  any  curve  starting  in  region  P 
will  have  -  0  while  any  point  curve  starting  in  region  E 
will  have  a  P^  -  1.  Por  regions  A,  B,  C,  and  D,  If  a  curve 
starts  within  ttie  limit  cycle.  It  will  spiral  towards  the  origin, 
hence,  P^  -  1.  If  It  la  outside  the  limit  cycle,  the  integral 
curve  will  diverge  from  the  limit  cycle  and  the  value  of  P^ 
will  become  1  or  0  depending  on  whether  the  integral  curve  enters 
region  E  or  P  first. 

Figure  3  Indicates  the  isoprobablllstlc  lines  for  the  case 
of  W-1,  X«— 1.  Figure  4  indicates  the  Isoprobablllstlc 

lines  for  the  case  k  -  0,  X  ■  -  1. 
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4.  Computer  Technlquea 

Tlie  results  Illustrated  in  Section  3,  and  additional  results 
were  obtained  on  the  Datatron  Computer  at  Purdue  University, 

Since  r(t)  »  +  k  with  equal  prolablllty ,  the  Integral 
resolved  Itjelf  Into  a  sum. 


^n+l^^n+1 


n+1 


u,v) 


P  (x" 
n  n 


ynJU. 


yn*  *{-  "  '‘n  *  ><} ' 


X  was  chosen  as  +1  or  —  1  on  different  runs,  k  was  chosen 
as  1,  .1,  0  on  different  runs.  (a,h)  was  chosen  as  (l,l) 
and  (4,4). 

A  grid  of  2^^  01  points  (“^il  x  51)  was  selected  for  ranges 
of  u,v; 


(a,t)  -  (1,1), 

-  2.5  <  u  <  2.^ 

-  2.5  <  V  <  2.5, 
Au  -  Ziv  =  .  1 , 


( a , l )  »  (4,4), 
0  <  u  <  5, 

0  <  v  <  5, 

/)u  -  Av  -  .  1 , 


At  -  .1,  At  -  .1. 

Initial  values  of  were  selected  for  all  points  of  the 
grid;  Pq  -  1  for  all  nolnts  (x,y)  satisfying  x  <  a,  y  <  b; 
Pq  -  0  for  all  other  points. 
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The  argximents  for  equation  (5')  were  computed  In  a  standard 
fashion,  and  probabilities  assigned  by  linear  Interpolation  over 
the  four  adjacent  comers  for  points  within  the  grid,  and  by  a 
like  extrapolation  for  points  falling  outside  the  grid. 

N 

^n+1  were  based  upon  a  full  grid  of  values. 

When  all  2601  values  were  calculated,  the  entire  grid 

function  was  then  replaced.  This  procedure  constituted  one 
cycle . 

On  the  Datatron,  one  cycle  required  an  average  of  70 
minutes . 

Thirty  cycles  were  calculated  for  the  Illustrative  figures 
^  above. 


BLANK  PAGE 


